-- ********************************************************************
--    $Id: 0.2 upgis_cascadeunion_functions 2008-10-05 $
--   ********************************************************************
-- 	 Copyright 2008, Regina Obe, Paragon Corporation
--   License: BSD
-- 	 Purpose: This is a set of plpgsql/sql functions  to implement aggregate spatial union
-- 	 It is a drop in replacement for the built in PostGIS ST_Union aggregate function but in general 5-10 times faster
--	 It is loosely based on JTS 1.9 cascade union implementation which will hopefully be available in the next version of GEOS
-- Usage Example:  SELECT somefield, upgis_cascadeunion(the_geom) As uniongeom FROM sometable GROUP BY somefield
-- ************************************************************************************ 
CREATE OR REPLACE FUNCTION upgis_dump_collect_garray(geometry[])
  RETURNS geometry AS
$$
	SELECT CASE WHEN $1 IS NULL THEN null::geometry
		WHEN array_upper($1,1) = 1 THEN $1[1] 
		WHEN array_upper($1,1) > 1 THEN st_collect_garray(ARRAY(SELECT CASE WHEN ST_NumGeometries($1[i]) > 0 THEN (ST_Dump($1[i])).geom ELSE $1[i] END FROM generate_series(1, array_upper($1,1)) As i))
		ELSE null END;
$$
  LANGUAGE 'sql' IMMUTABLE;

CREATE OR REPLACE FUNCTION upgis_unitecascade_garray(ginitset geometry[], iter int,
startcount int, nonintercompleted boolean)
  RETURNS geometry AS
$$
DECLARE 
	gcount integer := array_upper(ginitset,1); --initial count
	result geometry := null; -- final result
	gnonint integer[] := null; -- non-intersecting geoms index
	gstep integer := CAST(ln($3 + 1)/(ln($2 + 1)) As int); --max number of geometries to unite in any iteration
BEGIN 
	IF gcount < 2 THEN --If less than 2 geometries - return it no need to union any further
		result :=  ginitset[1] ;
	ELSIF gcount < 5 OR gcount < gstep  --If fewer than 5 or lower than step increment - union 
			or gstep < 2   THEN 
		result := st_collect_unite_garray(ginitset);
	ELSIF NOT nonintercompleted AND gcount < 1000 THEN
	--If fewer than 1000 records to look at and non-intersection check never run
	-- break out non-intersecting regions from intersecting and then collect non-intersecting
	-- with union of intersecting regions

	--Step 1: Determine index of all geometries whose bbox do not intersect any other
		gnonint := ARRAY(SELECT s 
			FROM generate_series(1, gcount) As s
			WHERE ginitset[s] IS NOT NULL 
				AND NOT EXISTS(SELECT ginitset[s2] 
					FROM generate_series(1, gcount) AS s2
					WHERE ginitset[s2] IS NOT NULL AND s <> s2 AND ginitset[s] && ginitset[s2])
			);
		--If count of non-intersect is 0 or null - then all intersect with at least one other geom cascade unite all
		IF gnonint IS NULL OR array_upper(gnonint, 1) = 0 THEN
			result := upgis_unitecascade_garray(ginitset, iter + 1, startcount, true);
		ELSIF array_upper(gnonint,1) = gcount THEN -- none intersect so just collect
			result := upgis_dump_collect_garray(ginitset);
		ELSE
		--some intersect and some don't intersect - leave all that don't intersect as separate geoms
		--cascade unite all that do intersect - and then collect the group
			result := upgis_dump_collect_garray(
				st_geom_accum(ARRAY(SELECT ginitset[gnonint[s]] As geom
				FROM generate_series(1, array_upper(gnonint,1)) As s
					
				), 
				upgis_unitecascade_garray(ARRAY(SELECT ginitset[s] As geom
					FROM generate_series(1, gcount) As s
						WHERE NOT (s = ANY(gnonint))
						), iter + 1, startcount, true)
				)
			); 
		END IF;
	ELSE
		
	--Otherwise just recursively union
		result := upgis_unitecascade_garray(ARRAY(SELECT st_collect_unite_garray(ginitset[i:least(i
			+ gstep - 1,gcount)]) As geom
			FROM generate_series(1, gcount,gstep) As i), iter + 1, startcount, false) ;
	END IF;
	return result;
END;

$$
  LANGUAGE 'plpgsql' IMMUTABLE;




CREATE OR REPLACE FUNCTION st_collect_unite_garray(ginitset geometry[])
RETURNS geometry AS
$$
DECLARE 
	gcount integer := array_upper(ginitset,1); --initial count
	result geometry := null; -- final result
	gint integer[] := null; -- intersecting geoms index
BEGIN
	IF ginitset IS NULL THEN
		result := null;
	ELSIF gcount < 2 THEN
		result := ginitset[1];
	ELSE
		--Determine index of all geometries whose bbox intersects at lease one other
		gint := ARRAY(SELECT s 
			FROM generate_series(1, gcount) As s
			WHERE ginitset[s] IS NOT NULL 
				AND EXISTS(SELECT ginitset[s2] 
					FROM generate_series(1, gcount) AS s2
					WHERE ginitset[s2] IS NOT NULL AND s <> s2 AND ginitset[s] && ginitset[s2])
			);

		--If count of intersect is same as count of whole set then they must all intersect at least once - just unite
		IF array_upper(gint, 1) = gcount THEN
			result := st_unite_garray(ginitset);
		--None of the geometries in this set intersect so we can just collect
		ELSIF gint IS NULL OR array_upper(gint,1) = 0 THEN --none intersect so we can just collect
			result := upgis_dump_collect_garray(ginitset);
		ELSE
		--some intersect and some don't intersect - leave all that don't intersect as separate geoms
		--unite all that do intersect - and then collect the group
			result := upgis_dump_collect_garray(
				st_geom_accum(ARRAY(SELECT ginitset[s] As geom
				FROM generate_series(1, gcount) As s
					WHERE NOT (s = ANY(gint))
				), 
				st_unite_garray(ARRAY(SELECT ginitset[s] As geom
					FROM generate_series(1, gcount) As s
						WHERE s = ANY(gint)
						)
					) 
				)
			); 

		END IF;
	END IF;
	RETURN result::geometry;
END;
$$
	LANGUAGE 'plpgsql' IMMUTABLE;


  
CREATE OR REPLACE FUNCTION upgis_unitecascade_garray_sort(geometry[])
  RETURNS geometry AS
$$
	SELECT CASE WHEN array_upper($1,1) < 5  THEN
st_unite_garray($1) ELSE
		upgis_unitecascade_garray(ARRAY(SELECT $1[i] FROM generate_series(1,array_upper($1,1)) i  WHERE $1[i] IS NOT NULL ORDER BY $1[i]) , 1, array_upper($1,1), false) END;
$$
  LANGUAGE 'sql' IMMUTABLE;


  
DROP AGGREGATE upgis_cascadeunion(geometry);

CREATE AGGREGATE upgis_cascadeunion(geometry) (
  SFUNC=st_geom_accum,
  STYPE=geometry[],
  FINALFUNC=upgis_unitecascade_garray_sort
);





